% Copyright (c) 2012, John Gebbie. All rights reserved.
%
% This library is free software; you can redistribute it and/or modify it
% under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation; either version 2.1 of the License, or (at
% your option) any later version.
%
% This library is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
% General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this library; if not, write to the Free Software Foundation,
% Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301  USA


% Script for testing convergence zones.


%% Init

clear;
clear classes;
clc;

obj = ortrace();

%% Parameters

max_range = 100000;
seabed_depth = 6000;

params = [];

% SSP

% Munk SSP
munk_z = (0:10:seabed_depth).';
munk_epsilon = 0.00737;
munk_ztilde = 2*(munk_z-1300)/1300;
munk_c = 1500*(1+munk_epsilon*(munk_ztilde-1+exp(-munk_ztilde)));
munk_r = 0;

iRng = 1;
params.ssp_cz{iRng}.c = munk_c;
params.ssp_cz{iRng}.z = munk_z;
params.ssp_r(iRng) = munk_r;
params.ssp_step_z = 100;

% SEABED
params.seabed_range = 0;
params.seabed_depth = seabed_depth;
params.seabed_c = 1600;
params.seabed_rho = 1.6;

% RAY
params.rs = 0;
params.zs = 30;
% ray angle (positive angles points toward seabed)
params.theta_mn_mx_n = [ [1 1].*0*pi/180, 1 ];
params.output_ranges_mn_mx_n = [0 max_range 100];
params.output_full_ray = false;
params.output_grazing_angle = false;
params.step_m = 25;
params.max_db_loss = 200;
params.max_depth = Inf;

%% Configure

tic
obj.init_environment(params);
toc

%% Cast Ray

tic
obj.cast_ray_fan();
toc

%% Plot

ray_out = obj.rays_out{1};

clf;

subplot(1, 4, 1);
xx = obj.ssp_c;
yy = obj.ssp_z;
plot(xx, yy);
set(gca, 'ydir', 'reverse');
axis tight;
ylim([0 max(params.seabed_depth)]);
title('c');

subplot(1, 4, [2 3 4]);
xx = ray_out(ortrace.POS_R, :);
yy = ray_out(ortrace.POS_Z, :);
plot(xx, yy, '-ob');
set(gca, 'YDir', 'reverse');
set(gca, 'YTick', []);
xlim([0 max_range]);
ylim([0 max(params.seabed_depth)]);
title('ray');

hold on; % draw seabed and sea surface
plot([0, max_range], [0 0], 'k');
xx = [ 0 params.seabed_range(1:end) max_range ];
yy = params.seabed_depth([1 1:end end]);
plot(xx, yy, 'k');
